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FOREWORD 


The Mul tispectral Resource Sampler (MRS) “Proof-of-Concept" 

Study is intended to be a comprehensive analysis of the corrections that 
must be applied to MRS data to allow for atmospheric correction factors 
and the variability of bidirectional reflectance from the scene. 

This study was initiated by Dr. Charles Schnetzler of NASA Goddard 
Space Flight Center, and was performed by ORI, Inc. Space Data and Systems 
Division with Mr. Charles W. Aitken coordinating the efforts of ORI 1 s 
consultants. 

The complete study results are reported in five separate volumes 
which have the following titles and authors: 

DETERMINATION OF ATMOSPHERIC OPTICAL PARAMETERS USING 
THE MULTISPECTRAL RESOURCE SAMPLER, by Dr. Robert E. Turner, 
Science Applications, Inc. 

ATMOSPHERIC CORRECTION USING AN ORBITAL POINTABLE IMAGING SYSTEM 
By Dr., Philip N. Slater, University of Arizona 

A SIMULATION ANALYSIS OF BIDIRECTIONAL REFLECTANCE PROPERTIES 
AND THEIR EFFECTS ON SCENE RADIANCE— IMPLICATIONS FOR THE 
MRS, By Dr. James A. Smith, Colorado State University 

'MRS LITERATURE SURVEY OF BIDIRECTIONAL REFLECTANCE, By 
Dr. James A. Smith, Colorado State University 

MRS LITERATURE SURVEY OF ATMOSPHERIC CORRECTIONS 
By, Dr. Philip N. Slater, University of Arizona 
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1.0 INTRODUCTION 


Methods for the correction of Landsat data for astmopheric transmittance 
and path radiance have thus far relied on the use of atmospheric models and/or 
extensive ground-based measurements. The use of an orbiting pointable imager 
may allow atmospheric corrections to be made without any supporting ground 
measurements or any knowledge of the ground reflectance. If a pointable 
imager were launched with the Thematic Mapper (TM) and/or the Multispectral 
Scanner System (MSS), or if it were flown in close orbital proximity to such 
sensors, then it might be possible to correct the data, collected by these 
systems for atmospheric effects for any area of the earth's surface. 

The following preliminary study was conducted to consider the feasibility 
of such a correction method because of its potential importance in improving 
scene classification accuracy, correcting bi-directional reflectance distri- 
bution function (BRDF) data, correcting obligue angle data collected for bio- 
mass estimation in sparsely vegetated areas, and for use in the analysis of 
the results of polarization experiments. 
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PROJECT OUTLINE 


1.1 


The work performed as part of this project has included: 

1) A literature survey of work in detemring the influence of atmospheric ef- 
fects on spectal signatures and classification accuracies. Over sixty 
references were listed and six chat were particularly relevant to classi- 
fication accuracy were discussed in detail. This survey was described in 
a report entitled "MRS Literature Survey of Atmospheric Corrections". 

(MRS is an abbreviation for Multispectral Resource Sampler), 

2) The spherical geometry of earth sensing from a pointable sensor in a 
Landsat orbit was examined. This is of value for work in both bi- 
directional reflectance distribution function (BRDF) studies as well 

as atmospheric correction. The ranges of solar aximuthal zenith angles, 
and Air mass which can be encountered as a function of latitude and time 
of year are described as a function of pointing angle. 

3) The literature survey, mentioned in item (1) above, revealed that changes 
in atmospheric extinction coeffecient (or optical thickness), x ex ^, 
greater than 0.1, between training and test site, began to have a signifi- 
cant effect on classification accuracy. A study is reported here that was 
made to corroborate this result and determine the sensitivity of classifi- 
cation accuracy to change in atmospheric effects. 

4) A two-wavelength, single target method for determining t was investigated 
in detail but was found to underestimate T gxt because the rate of change 
in path radiance (L p ) with pointing angle, 6 , was greater than the rate 
of change in total radiance (L ) at the sensor, 

5) A non-linear regression method is described which enables T fiXt , path 
radiance, Lp, the reflectance, p, of a uniform Lambertian ground feature 
to be determined from a knowledge of the variation of the radiance at 
the sensor, L & , with pointing angle e, for a particular wavelength 
interval AX. 
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Plot of solar azimuth angle veicsus solar declination angle 
for a typical Landsat orbit. ‘Flie day and month of die year 
and the solar zenith angle are marked for 5° increments of 
declination angle. The latitudes considered are 35°W, 40°N 
45°N and 55°N. 






The size of the Instantaneous Field of View (IFOV) projected on to the 
ground as a function of pointing angle is important because it indicates the 
size of the uniform ground reference area required for accurate atmospheric and 
BRDF measurements. Figure 3 shows the increase in size of an IFOV 15 m X 15 m 
at the ground at nadir as the pointing angle to nadir increases. The easy con- 
sidered Is that for an orbital altitude of 700 km and the effect of earcb cur- 
vature has been taken into account. It is worth noting for comparison purposes 
that the IFOV length equals that o* the MSS IFOV (76 m) at an angle of about 
56°. The pointing direction becomes tangential to the earth's surface at a 
pointing angle to nadir of about 66°, 

The position of the sensor in its orbit with respect to a reference area 
on the ground for various pointing angles is shown in Figure 4. Again the or- 
bital altitude is 700 km. The adjacent table indicates that a maximum air mass 
of 4 is encountered where the pointing angle is nearly tangential to the earth's 
surface, at 65°, in orbital position E. 

The general atmospheric problem in remote sensing is illustrated in Fig- 
ure 5. The figure shows radiant flux from the sun incident on two ground areas 
B and A and a remote sensing system pointed at area B. As the direct flux from 
"un passes through the atmosphere to the ground it is attenuated by scatter- 
ing and absorption. (We shall assume here that we are sensing through non-ab- 
sorbing spectral regions of the atmosphere.) The ground area B is irradiated 
primarily by direct solar flux partially attenuated by scattering such .as in- 
dicated at C. In addition there is a diffuse component of the irradiance due 
to scattering in the atmosphere shown schematically as originating from G but 
actually coming from the entire hemisphere over the ground area. Flux reflected 
from B is attenuated by scattering, such as at F, on its path from the ground to 
the remote sensing system. In addition scattered direct solar flux, as at D, is 
directed into the entrance pupil of the sensor. There is also a crosstalk compo- 
nent owing to flux reflected from area A being scattered in the direction into the 
sensor entrance pupil as at E. These last two scattered flux components are 
indistinguishable from the flux coming from the surface of interest and are com- 
monly referred to as path radiance. (Owing to multiple scattering effects in 
the atmosphere there are other components of path radiance, however the two 
shown in Figure 5 are in general the largest contributors) . A point wo.: th noting 
with respect to Figure 5 is that the magnitude of the components contributing 
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Figure A. A diagram showing the position of a pointable sensor 
with respect to a fixed ground scene as a function 
of pointing angle. The adjacent table shows the 
increase of air mass with increase in pointing angle, 
(Slates, *980). 
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to path radiance change with pointing angle. For example, as the pointing 
angle is increased from nadir to near-tangential , the radiance of the ground 
surrounding the reference area will provide an increased contribution to the 
path radiance. 

The desirability of a large change in air mass between the ground and 
sensor in order to have a large variation in sensor radiance to measure is 
tempered by the above considerations as well as a number of engineering 
factors. These are summarized below: 

1. If horizontal inhomogenieties in atmospheric conditions exist in the 
region of the reference area, the interpretation of the sensed data 
will become more difficult as the pointing angle increases. 

2. The greater the pointing angle, the larger the ground-projected IFOV 
and the larger the uniform ground reference area has to be. 

3. The influence of the radiance of the surround on the total radiance at 
the sensor increases with increase in pointing angle. This "crosstalk'’ 
effect will complicate the interpreation of the data and reduce the 
accuracy of the measurement results. 

4. Large pointing angles require a larger and heavier pointing mirror with 
a greater moment of inertia than otherwise necessary, 

5. The larger the continuous pointing angle sweep made per target, the 
fewer the number of observations possible. Reference to Figure 6 shows 
that a sweep in pointing angle from 63° to 0° takes about 300s, which is 
about the time for the spacecraft to cross the forty-eight states. In 
otherwords, if such a large continuous angular sweep is required, no 
other near-nadir observations could be made of the United States during 
that particular orbital pass. In comparison, a 45° pointing angle 
sweep takes 100s. 
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3.0 SENSITIVITY OF CLASSIFICATION ACCURACY TO ATMOSPHERIC CHANGES 


A conclusion from the literature survey conducted at the start of this 
project and reported in "MRS Literature Survey of Atmospheric Corrections" 
was that changes in the atmospheric extinction coefficient x ext greater 
than 0.1 causeu a substantial decrease in classification accuracy. The work 
described in this section was carried out to independently verify these 
earlier results using a different approach and to develop a more complete 

understanding of the relationship between changes in x ext and scene clas- 

sification accuracy. 

To gain an appreciation for the magnitude of x gxt at different sites 
across the US, we refer to the work of Flowers et al . (1969) who determined 
the cumulative frequencies of the daily average turbidity for typical 
urban, suburban and rural conditions. Note that the turbidity coefficient, 
B, used by Flowers et al. is the decadic extinction coefficient at a wave- 
length of 0,5 ym and is related to the often-used atmospheric extinction 
coefficient, or optical depth, x, by B = 2.3x. x is often quoted for X = 

0.55 ym whereas Flowers, et al . , used B values at 0.5 ym. The effect of 

the change in wavelength depends on Rayleigh scattering and the atmospheric 
aerosol distribution; however, as a first approximation Xq 55 = Q.9xg 5 . 

The plot of Flowers, et al . , for sites across the United States, of the 
cumulative frequencies of the daily average turbidity for typical urban, 
suburban and rural conditions is presented in tabular form in Table 1 in 
terms of optical depth x for a wavelength of 0.5 ym. (The subscript 'ext' 
has been omit.ted). Note that 95% of the days in rural areas huve x values 
of less than 0.3 which corresponds in a standard atmosphere, to a visibility 
or meteorol igical range of 25 km. 

If only changes in x of greater than 0.1 across a scene significantly 
reduce classification accuracy, then the results of Flowers, et al . , show 
that the classification of purely rural areas will usually be independent 
of atmospheric conditions. However, in other cases, for example in land use 
studies at urban - rural interfaces, changes in x across the scene can be 
expected to be much larger than 0.1. Thus the table of Flowers, et al., 
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TABLE 1 • 


PERCENT OF 

DAYS 

WITH HAZE 

LEVEL IN 

INDICATED 

RANGE 



Lower/upper bounds for 

T 


0/0.1 

0. 1/0.2 

0.2/0. 3 

0. 3/0.4 

0 . 4/on 

Rural 

10 

65 

.20 

4 

1 

Suburban 

7 

28 

35 

12 

18 

Urban 

2 

13 

15 

15 

55 


shows that 65% of the days in a typical rural area have t values between 
0.1 and 0.2 and 55% of the days in a typical urban area have x values 
greater than 0.4. 

In this study of the influence of the change in atmospheric conditions 
on classification accuracy and later, in the study of methods to determine 
atmospheric parameters from measurements made by a pointable imager, atmos- 
pheric data were derived from a radiative transfer calculation technique 
described by Herman and Browning (1975). This technique utilized model 
atmospheres constructed by Elterman with various amounts of aerosol loading. 
The model assumed a flat earth, a uniform Lambertian surface and a hori- 
zontally homogeneous atmosphere. All orders of multiple scattering were 
considered in the computations. Polarization data are available from this 
calculation but were not used in these preliminary investigations. The out- 
put from the radiative transfer computation gave the total radiance and path 
radiance emerging from the top of the atmosphere and the total irradiance at 
the ground for various viewing azimuth angles and view angles to nadir (point 
ing angles) as a function of ground reflectance, solar zenith angle, wave- 
length and atmospheric extinction coefficient, 

A range of the atmospheric spectral extinction coefficients used as 

input to the Herman and Browning model are shown in Fig. 7. A Rayleigh 
atmosphere is shown as the lowest curve, the next is a very clear atmos- 
phere typical of the desert southwest on a day of greater than 100 km 
visibility, The next three curves are for x values at 0.65 ym that increase 
in steps of 0.05. The top curve is for very poor meteorol igical conditions, 
a x value of 0.65 at 0.5 ym. It is unlikely that classification work 
would be attempted under atmospheric conditions as poor as this. 

The remote sensing data used to provide the classification for the 
generation of simulated data were Landsat MSS Computed Compatible Tape 
(CCT) values taken from an image of the San Carlos Indian reservation, an 
arid, rural region in central Arizona. The meteorological conditions were 
excellent at the time the imagery was collected, corresponding to the curve 
labeled ‘very clear 1 in Fig. 7. (This is an estimate as no atmospheric 
measurements were taken at the time the imagery was collected and only 
qualitative r< oorts of the atmospheric condition are available). Simulated 
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data were generated from the Landsat CCT in four steps: 

1, Three surface classes were selected from the imagery: limestone, 
riparian vegetation and basalt. These were selected because they cover 
a reasonably comprehensive range of reflectances and spectral character- 
istics. The riparian vegetation was fairly dense (cover > 60%) but 

the relatively high reflectance values in bands 4 and 5 are indicative 
of the influence of underlying soil and rock on the signature. The 
means and standard deviations of the gray levels for the limestone, 
vegetation and basalt areas are shown in Fig. 8 for the four MSS bands. 
(Note that the band 7 values were multiplied by 2 to account for the 
difference in the quantization range of the data (0 - 63 instead of 
0 - 127). 

2. The simulated data were first generated as normally distributed inde- 
pendent random numbers with zero mean and unit variance, independently 
in each of the four spectral bands for each class. 

3. Using the eigenvector and eigenvalue matrices and mean vectors of the 
specified classes, these data were transformed to new data with the- 
desired mean and covariance matrices of the Landsat test classes. 

Since all the transformations were linear, the final data were also' 
normal independent random, variables. Table 2 lists mean, variance, and 
covariance and correlation matrices for the basalt, limestone and vege- 
tation training classes. Each class consisted of 2400 pixels, quantized 
to 128 gray levels. 

4, Using path radiance values generated from the Herman and Browning at- 
mospheric model, the synthetic image data were transformed from the 
originally assumed x value of 0.125 at 0.55 ym to t values of 0.25 and 
0.56. An important point to note here is that the difference between 
the original and the transformed data is purely deterministic . Thus 
there is no inclusion of natural deviations between training and test 
sites. The classification analysis depends only on the different 
atmospheric conditions assumed. 

Fig. 9 is a gray level picture of the synthetic data set for band 4. 

The top three pictures are for the very clear atmosphere, t = 0.125, and 

the bottom three pictures are for the t s 0.56 case. Note how the contrast 
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Figure 8. The mean and standard deviation gray 
levels for limestone, vegetation and 
basalt as recorded by the four Landsat 
MSS bands under very clear atmospheric 
conditions. 
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COVARIANCES, MEANS, COVARIANCE MATRICES AND CORRELATION MATRICES FOR HIE THREE GROUND FEATURES. 
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Figure 9. A band 4 gray level picture of a 

typical data set. The top is for a 
very clear atmosphere and the bottom 
is for a very hazy atmosphere. 
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of the pictures decrease as the atmosphere becomes more hazy. This effect 
is less noticeable in the case of limestone because it has a higher re- 
flectance than the other surfaces. 

The class statistics for the three atmospheric conditions are shown in 
Table 3. Note that the class means increase with increasing haze and that 
the increase is most noticeable for the lower reflectance surfaces. The 
standard deviations decrease as the atmosphere becomes more turbid. This 
again is due to a reduction in contrast or a compression in the gray level 
range. 

LARSYS was used to classify the atmospherically degraded imagery. 

Figure 10 shows the result of training on the features under very clear 
atmospheric conditions and testing under hazy conditions. The errors in 
classification are seen to be greatest for basalt, This is due to basalt 
having the -lowest reflectance of the three features and therefore its sig- 
nature being modified most by changes in atmospheric conditions. 

A second classification was made using LARSYS in which the training was 
performed under the hazy atmospheric conditions and the tests were conducted 
under the clear atmospheric conditions. The results are shown in Fig. 11 
and it can be seen that the errors in classification, for the same change 
in atmospheric conditions as in the previous case, are less than those in 
Fig, 10. 

The pictorial results illustrrted in Figs. 10 and 11 are summarized 
quantitatively in Table 4, along with the results of other atmospheric 
changes. The first numerical column indicates the accuracy of the training 
site area. Three important conclusions can be drawn from Table 4: 

1. The decrease in classification accuracy is greater the greater the 
change in atmospheric conditions unless the initial accuracy is 
very high when the accuracy can actually increase (see the limestone 
entry in Table 4 ). 

2. The classification accuracy decreases the most for basalt, the lowest 
reflectance feature. 

3. Except in the case of limestone, the classification accuracy de- 
creases less with change in atmospheric conditions when the training 
is performed under poor atmospheric conditions and the testing is 
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Band Atm 0 


Basalt 

4 31. 9 

5 48.4 

6 52.1 

7 3 


Vegetation 


4 30.7 

5 41 .0 

6 63.8 

7 60. 1 


Limestone 

4 53.1 

5 77.2 

6 82.4 

7 67.1 




}! 


TABLE 3 

CLASS STATISTICS 


Mean Standard deviation 


Atm 1 

Atm 2 

Atm 0 

Atm 1 

Atm 2 

33.8 

37.4 

2.46 

2.41 

2.22 

49,1 

51.2 

4.73 

4.81 

4.65 

52.9 

52.9 

5.21 

5.15 

5.07 

41.9 

42.5 

5.41 

5.44 

5.21 


32.6 

36.2 

5.61 

5.39 

5.07 

41.8 

44.2 

11.59 

11.64 

11.01 

64.4 

64.4 

4,97 

5.03 

4.96 

60.7 

60.4 

6.58 

6,59 

6.24 


53.7 

56.0 

5.65 

5.61 

5.13 

78.0 

79.0 

8.42 

8.43 

7.86 

83.2 

82.0 

9.06 

8.99 

8.62 

67.6 

67.0 

7.79 

7.81 

7.37 
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Figure 10. 'The result of a LAilSYS classification using a clear atompshere (top) for 
the training conditions and a hasy atmosphere for the test conditions. 
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figure 11 . The result of a LAR 3 Y 5 classification using a hazy atmosphere ( top ) for the 
training conditions and a clear atmosphere for the test conditions 
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CLASSIFICATION ACCURACIES 



ATMP/ATMP 

ATM0/ATM1 

ATMJ3/ATM2 

ATM1/ATM2 

BASALT 

98 , it 

93.5 

44.7 

91 .8 

VEGETATION 

96.1 

95-7 

89.3 

92.0 

LIMESTONE 

98.6 

99.5 

100.0 

99.8 

AVERAGE 

97.7 

96.2 

78.0 

94.5 


ATM2/ATMI 

ATM2/ATMJ3 

ATM1 /ATMP 


BASALT 

98,4 

97.8 

99.5 


VEGETATION 

95.0 

93.3 

94.3 


LIMESTONE 

93.8 

92.7 

97.2 


AVERAGE 

95.7 

94.6 

97.0 
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under good atmospheric conditions than vice vem. This phenomena has 
been reported by Potter and Shelton (1974) and Turner (1975). 


26 


ORIGINAL PAGE IS 
OF POOR QUALITY 


4.0 THE TWO-WAVELENGTH APPROACH TO ATMOSPHERIC CORRECTION 


In order to correct the spectral radiance at the entrance pupil of a 
sensor for atmospheric effects, we need to know the transmittance along the 
path from the ground feature to the sensor and the path radiance, l.e,, 
the radiance detected by the sensor introduced by atmospheric scattering. 

In this and the next section we shall describe two different approaches to 
determining the atmospheric extinction coefficient, x (related to 
transmittance, x, by x = exp - x ext ) and path radiance Ip. 

The basis of the two-wavelength approach is as follows. If tha total 
6 dependent radiances at the sensor at, for example, wavelengths 0.69 and 


Ly-| , the ground scene is Lambertian with a change in 


0.71 urn are L gg and 
reflectance from p at 0.69 ym to p + Ap at 0.71 ym then 


■69 


= Eo 

7T 


=x sec0 


. Ax sec©\ 
e ; 


+ L p (e) 


and L 


71 


E(p + Ap ) 

7T 


^X 


sec© e ”Ax sec0l 


+ Lp(©) 


where t is the value of the extinction coefficient at 0.70 ym. The ground 

irradiance at the neighboring wavelengths is assumed the same. We can then 

subtract these equations to obtain \ 

i . _ Ep rpr Ax sec -Ax sec© - Ap „ -Ax sec©) 

len - L „ — - e e -e — - e 


■69 


71 


Taking the natural logarithnrof both sides 


ln(L 6g - L^ = -x sec© + In + In ( ) 

where the last brackets include the three terms in parentheses in the pre- 
vious equation. Obviously, if the last two terms in this equation are 
independent of sec© then the slope of a plot of In (L gg - L^) versus sec© 
will yield 7 . 

A typical result of this method is shown in the lower curve in Fig. 12. 
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Eigure 12. A plot of In (L q ^ - L Q 
versus secQ, 

upper curve corrected for path 
radiance lower curve uncorrected 
for path radiance. 
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Model data were used of x * 0.110 at X * 0.69 pm and t b 0.105 at X * 

0.71 ym. A change in reflectance, representative of vegetation, of 0.05 
from 0.05 to 0.10 was also chosen. The curve of In (Lgg - L^) versus sec© 

is a straight line with all the points falling close to the curve. Unfor- 
tunately, the slope of the curve yields a low x value of 0.062. In order 
to determine the cause of the discrepancy in x value, th^ value of Lp 
calculated from the model output were subtracted from the total radi- 
ance value of l p at the sensor and the above procedure was repeated. The 
upper curve in Fig. 12 resulted, which gave a x value very close to the 
correct mean value. This indicated that the problem lay in assuming the L p 

values increased at the same rate with 0 as the total radiance values at the 

# 

sensor. If in fact the L p contributions did not cancel with subtraction, 
the total radiance and the path radiance for the two wavelengths and re- 
flectances as a function of © are shown in Figure 13 and the difference in 
the slopes between the L s and Lp curves is evident. Unfortunately, although 
the choice of reflectance and extinction values yielded a difference in -L s 
values twenty times larger than the difference in Lp values, the greater 
slope of Lp versus sece strongly influences the estimated value for x .’ 

Other values of reflectance at X = 0.69 and 0.71 p m were tried and the 
results are listed in Table 5. These show that the values of the reflectances 
used do affect the accuracy of the estimate of x but this influence is over- 
shadowed by the ’slope of the curve of L p versus sec6. 


L RELATIVE TO AN E 0 VALUE OF 1000 



t VALUES 
T 690 

P 710 


TABLE 5 


0.110 


* 71 ° ' ° 

.105 



P 690 



0.05 

0.1 

0.15 

0.1 

0.15 

0.20 

0.061 

0.050 

0.046 

0.058 

0.049 

0.057 


t values 
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5.0 A NUMERICAL PROCEDURE FOR DETERMINING x, p, and l 

r 

The acquisition of radiometric measurements via a pointabl e-satellite- 
borne radiometer would provide nadir-angle dependent data, and consequently 
more information about the atmosphere. A numerical procedure to analyze 
these measurements will be discussed that will permit determination of optical 
depth (x), path radiance (L p ), and ground reflectance, (p). This method 
simultaneously computes these parameters, effectively applying the atmos- 
pheric corrections during processing. 

Mathematical Formulation 

A non-linear regression procedure was utilized to fit nadir-angle de- 
pendent radiometric data to a theoretical expression describing the measure- 
ments. Radiative transfer calculations provided realistic data to test this 
method. 

In general, the radiance at the entrance aperture of a pointable satel- 
lite-based radiometer for a specified spectral band is given by the ex-* 
pression 

l s (s N ,e z ,f,p) = E tot (0 N ,a z ,T,p) ^yz 1 expfcsec (© N ) ) + L p (0 N> e z ,T,pj (1) 

where the functional dependencies of the quantities are indicated. 

The quantities involved are: 

L s : Radiance at the sensor 

E to t : Total irradiance at the ground 

E tot - E Q cos (0 Z ) exp (-xsec (0 Z )) + E Q (2) 

E q : Extraterrestial solar irradiance in spectral bands 

E d : Diffuse sky irradiance at ground 

x : Total atmospheric optical depth 

0^ : Nadir angle 

0 Z : Solar zenith angle 

p : Ground reflectance 

Lp : Path radiance 
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The quantities t, p, and Lp need to be known in order to analyze the 
ground scene. Unfortunately, this is a problem of at least four unknown 
Ep » x, p, and L p , with and Lp having an additional angular dependence. 

A series of measurements of as a function of 0^ will not analytically 
resolve this inadequately specified problem, except for the nonreallstic 
case of x = 0 where Lp - 0 and E D * 0. Combining information from various 
spectral bands may produce an analytic result with less severe constraints 
on t, but would require knowledge about the spectral character of Eg, t, p, 
and Lp. These spectral characteristics, however, could only be derived from 
information about or aprlorl assumptions regarding the atmosphere. 

Given a set of L g (©^) values, the unknowns can be determined by a non- 
linear regression technique if the angular dependencies of p and Lp can be 
specified. Model atmosphere radiative transfer calculations of L g and L p 
revealed a simple functional dependence for L p . The radiative transfer 
analysis (Herman, 1975) assumed a plane parallel atmosphere and a uniform 
Lambertian reflectance. For the representative conditions of a solar azi- 
muth angle of 90° and zenith angle of 45° the model calculations show L p to 
be a linear function of exp (-xsec.©^), Fig. 14. Therefore, for a given x, 

, p, and fixed © z , we can expand L , such that 

L p (© N ) = a - 6 exp (-xsec0 N ) (3) 

Note that this linear relationship holds only for solar zenith angles near 
90°. For other angles a simple polynomial expression (probably quadratic) 
is necessary. 

The angular variation of p could be included by parameterizing the 
appropriate bidirectional reflectance distribution function (BRDF) for the 
type of scene viewed. However, the ultimate goal in this analysis was to 
determine if the numerical technique could predict the radiative transfer 
calculation values for inputs x, p and outputs Eg, Lp . Since the radiative 
transfer computation assumed a constant p, a BRDF was not incorporated in 
the numerical analysis, but it could be accommodated. 

Therefore, with the above assumptions regarding the angular depen- 
dencies, equations 1-3 can be combined to give a general theoretical ex- 
pression and representation of the results of the radiative transfer 
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10 (REFERENCE UNITY E 0 , Q z = k5°, $ = 90°) 


p -.75 


i. 


•T sec0| 


Fig. 14. Plots of path radiance, L p , versus exp(-x sec6^) for 

X = 500 nm, T = 0.178, 0.643, and P = 0, 0.1, 0.25, 0.5 
0.75. Note the linear character of these plots for the 
choice of a solar zenith angle of 45° and azimuth angle 
of 90°. 
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calculations for the total radiance at the sensor as a function of nadir 
angle. 


MV 


a + 


* 

~ (E o cos© z exp(-T sec0 2 ) + Eg) - B exp(-t sec©^) 


( 4 ) 


5.1 Non-Linear Regression Technique 

Equation (4) can be rewritten as a non-linear expression of unknown 
parameters a.. for a range of nadir angles © 1 (the N subscript will be 
suppressed, but Implicit) 


MV ■ a 4 + 


\ (E q cus0 z exp(-a 1 sec0 z ) + a g expC-a^ecQ^) 


(5) 


where a 1 » t, a 2 » p, a 3 = E Q , and a 4 , a g determine L . Given a set of 

measured values, L^, as a function of ©.. and uncertainties, , in the 
data points L C1 - , a non-linear regression technique can optimize the choice 
of parameters a. by minimizing x*> where 

J 


2 r 
X 2 

i 


i [m - l s <vT 


f 

a i 


• ( 6 ) 


This estimate of the parameters, a^ , will determine the most probable 

functional form of L s (e^) that will fit the data L s ^ and consequently give 

an estimate for t, p, and L_. 

P 2 

The regression analysis estimates the parameters by minimizing x with 

respect to each of the parameters simultaneously, 



Since L e (e.) is a function of five parameters, this method finds the minimum 
of x in a five-dimensional hyperspace. The numerical technique utilized 
to accomplish this search of the hypersurface was first suggested by 
Marquardt (1963) and coded in Fortran by Bevington (1969). 


35 


The Marquardt approach probes the hypersurface by determining the 
gradient of the surface given an initial set of parameters, a^, and para- 
meter increments 5 a^. The a^ values are then Incremented following the 
path of steepest descent to the region of a minimum. Near the minimum the 
function L s (e.|) is expanded in a first-order Taylor's expansion, linear in 
the parameter increments fia^. Substituting the expansion into Eq. (7) 
and taking the derivative numerically generates a set of five simultaneous 
equations. The equations are then solved by linear regression matrix tech- 
niques to determine the optimum set of Increments 6a^ to augment the para- 
meter values to correspond to those at the hypersurface minimum. The 
process is repeated until the x value changes by less than a specified 
amount between iterations. The procedure then computes the uncertainties 
oaj, in the determined best fit parameter values, a^. 

Considering the above steps in detail, the first-order Taylor expansion 
in parameters a. can be written: 
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MV = + ? 


3L cn (e,) 

so' v< 


3a, 


(8) 


which is linear in 6 a • . L SQ and the derivatives are evaluated with the 
initial a., values. 


The derivatives are computed numerically for a given value 
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where Aa. is the non-optimized initial input guess for the parameter incre 

J 

ments 6a.. Substituting equation (8) into equation (6) yields 
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In order to minimize this linear regression problem with respect to the 
parameter increments 6a., expression (10) is differentiated, 
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such that for a given 6, , 

• Vi> • ic^r) “j- 
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Letting 3*. ■ - |^- * equation (11) can be recast as the matrix relation, 


where 
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jk 


B k * | Sa j “jk 
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Finally, solving the matrix relation 

g - 6a a 

for 6a will give the array of optimum increment values 3a. 


(14) 


(15) 


In order that this matrix method will quickly locate a minimum in the 

o 

X hypersurface, Marquardt suggested the modification 


S = 6a a' 


(16) 


where 


a' = a jk 
: a jk 


(1 + X) j = k 
0 t k 


(17) 


This alteration of the diagonal elements by the addition of X serves a 
useful purpose. For X small, 

a ~ a' 

fw *v 

•V ~ 

and the solution obtained for the . values is that of the Taylor 

J 

expansion. For X large, the diagonal terms dominate 


(18) 


6 . (19) 

J 

and the 6a. values increment a. along a scaled gradient toward the minimum. 
J j 
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Therefore, the non-linear regression algorithm combines a gradient search 
with an analytical solution developed from linearizing the fitting function 
given by 

6a = 0 e', (20) 


where c' is the inverse of matrix a 1 . 

*w A# 

A# 

Finally, the standard deviation of the computed parameters a- can be shown 

vJ 

to be related to the square root of the diagonal elements e' ... 

J J 

The algorithm is implemented in an iteration procedure. The program 

searches for a X value for the initial choice of a such that 

X 2 (a + 6a)<x 2 (6a). The array a is then set equal to a + 6a and the 

2 

analysis is repeated until X, or the change in x between iterations is 
less than a preselected value. 

The subroutine’ CURFIT (Bevington, p. 237) was utilized to perform a 
single X iteration. A listing of the FORTRAN computer code SKYRAD which 
utilized CURFIT in the regression technique analysis is contained in 
Appendix A. The program consists of the main program which input/output's 
the data, and calls the subroutines CURFIT, FUNCT, FDERIV, SCHISQ, and 
MXINV, and evaluates the results. 

The main program and subroutines are sufficiently documented to under- 
stand their operation. Comment cards in the main program describe the 
input data arrangement necessary to run the code. 

A brief description of the routines are: 

SKYRAD : Main Program 

CURFIT : Non-linear Regression subroutine described 

FUNCT : Evaluates Eq. (5) at for a given set of a^ values 

FDERIV : Evaluates the derivatives in Eq. (5) numerically via 

Eq. (9) at each point 0^ with respect to each parameter 
a.. 

J 2 

SCHISQ : Evaluates the reduced x value for a given set of data and 

fit parameters a. 

\) 

MXINV : Evaluates the inverse of the a 1 matrix. 

A typical output of the SKYRAD is contained in Appendix B. 
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Evaluation of the Effectiveness of the Regression Code 

In order to evaluate the efficiency and accuracy of the non-linear 
regression technique In computing the parameters before analyzing model 
atmosphere radiative transfer results, the program SVNRAD was written to 
synthesize data sets as Inputs to SKYRAD. For a given representative choice 
of a^ and 8^, SYNRAD computed the array L s (e.|) using Eq, (5) and the corre- 
sponding array of values for a specified amount of random noise added to 
L$ (9^ ) . The arrays 8^, L s (8), and cr^ were then stored as the input file 
TAPE5 for the code SKYRAD, 

The ability to exactly specify the input parameters of the synthe- 
sized data permitted examination of the following areas: 

1) Number and range of input 8^ values necessary 

2) Effect of input choices on the convergence and output of SKYRAD 

3) Independence of the parameters and their effect on L s (3.j} 

4) Effect of random noise on convergence 

5) Curvature of the x five-dimensional hypersurface. 

This test approach provided useful familiarity with the operation 
'of SKYRAD and the characteristics of the L (e) expression, Eq. (5). This 
understanding led to a procedure with which to analyze the L (e^ ) data from 
the model calculations. 

Graphical examination of Eq. (5) for various choices of the input 
parameters revealed that the p parameter strongly influences the shape and 
curvature of the L g (9^ ) data. This fact will be utilized in the next section. 
Four representative cases were consequently examined including the investiga- 
tion of the range of 8^ values. The parameter combinations considered were: 


X Uj) 


tp (a^, a g ) 

P (ag) 
0.1 

8,j Range 
0° - 60° 

0.2 

0.2 E 

o 

t— ■ 

cn 

m 

O 

0° - 45° 
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The synthetic data generated with these parameters were analyzed 
with SKYRAD for various choices of input parameters a.. 

J 

The choice of sixteen evenly spaced data points, 0^, gives the pro- 
gram eleven degrees of freedom with which to search for a solution that can be 
accepted with reaspnable confidence. This number of points also permits an 
adequate sampling of nadir values. There was no significant difference in 
SKYRAD outputs for identical cases of input parameters with either choice 
of 0.. range. Therefore, the restricted range, 0° £0^ £45°, was chosen be- 
cause for this range the ground-projected instantaneous field of view (IFOV) 
is smaller and the plane parallel atmosphere assumption of the atmospheric 
model radiative transfer computations is not seriously violated. 

To evaluate the effect of the first iteration initial guess of input 
parameters, a., a series of twenty-five input selections were evaluated for 

J 

each set of synthesized data. A range of values for a given a. were input 
holding the other parameters a k , jfk, constant at the modal values. Typically, 
each trial converged in three to ten iterations. 

The results of systematic parameter variation analysis are: 

a^: If t is overestimated, SKYRAD will converge to the correct 

solution but takes many iterations, 20-30. If x is under- 
estimated the code will converge, but a^ is changed to compen- 
sate and the predicted L p tends to be low. 

The solution SKYRAD predicts is very dependent on the choice 
of p. 

a 3 : The convergence of SKYRAD is very insensitive to variations in 
E D’ 

a 4 , a g : The SKYRAD output is very dependent on the input value of L p 

converging poorly for an overestimated L p input. 

a^: If all the parameters are systematically varied, either all 

substantially too large or too small, SKYRAD fails to converge 
or converges poorly. 
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Therefore, a judicious choice of input parameters is necessary to effectively 
use this numerical technique. 

The results of this investigation revealed significant details 

about the functional dependance of L.(0 4 ) on the a. parameters and the 
2 s i j 

curvature of the x hypersurface. The parameters of a^ are not strictly 

mathematically independent in their effects on the shape of L (6.). Vari- 

*9 I 

ations in one parameter can be compensated for by variations in other param- 
eters without significantly altering the curvature. In effect, the curvature 
of L vs 0. is not uniquely determined by a single choice of a. Consequently, 
the hypersurface is not characterized with a major deep minimum, but 
rather can be imagined as being bumpy with minor local minima, which SKYRAD 
can interpret as a solution. 

The fact that SKYRAD can converge to a number of equally acceptable 
solutions via the x 2 condition is disconcerting, but additional knowledge 

about the dependence of Ep and Lp on t and p enables a judicious evaluation 

of these solutions. The procedure necessary to effectively utilize the re- 
gression analysis will be discussed in the later section dealing with data 
from the model atmosphere radiative transfer calculation. 

An additional means of avoiding undesirable solutions was to prevent 
SKYRAD from predicting negative values for any of the a. parameters, since 

si 

only positive values are allowable given the form of Eq. (5). The subroutine 
CURFIT was modified to prevent determination of negative a. values for any 

J 

iteration. 

The effect of random noise added to the L s (e..) data was considered. 
The presence of up to 0.1% random noise did not prevent the convergence of 
SKYRAD to acceptable solutions. The presence of 156 random noise began to 
affect the fits to the synthesized data if po.0.5, when L (0^ ) data had little 
curvature. The net result of this analysis is to suggest that the instrument 
have a signal to noise ratio, -SNR > 100. 

Non-Linear Regression Analysis of Model Data 

In an effort to test the accuracy with which the numerical method could 
analyze real data, the output of radiative transfer calculations on model 
atmospheres was used as simulated "real" data input for the SKYRAD program. 


41 


SfflUSW 

Sets of (e.j, L s (6^)) data were compiled for specified wavelengths, but the t 
and p values were unknown to the processor. 

The radiative transfer computations were performed using the numerical 
approach described and coded by Herman and Browning. These calculations in- 
corporated various model atmospheres of varying amounts of aerosol loading 
as suggested by Elterman. 

For a fixed solar zenith angle of 45°, a unit solar irradiance, E , and 
a uniform Lambertian reflectance, the path radiance and the radiance at the 
sensor (L , L g ) were computed for a range of reflectances and nadir and 
azimuth angles. These calculations also permitted determination of Ep for 
each i, p combination. The computations were performed for wavelength values 
of 380, 500, 650, 750, and 900 nm. The optical depth for each wavelength was 
determined by the model atmosphere and amount of aerosol loading. A thin and 
thick aerosol optical depth were chosen for each wavelength. The major cases 
considered are listed below. 


X 380 500 650 750 900 

x RAY 

.450 

.145 

.026 

.0 

23 

.0 

L3 


x MIE 
x TOTAL 

.045 

.495 

.669 

1.119 

.033 

.178 

.498 

.643 

.048 

.074 

.397 

.445 

.027 

.050 

.344 

.371 

.020 

.033 

.300 

.313 



A typical computer output (X = 500 nm, x = 0.643) is contained in Appendix C. 

Figures 15, 16, and 17 are plots of the radiative transfer calculation 
outputs, L s (e.) is very reflectance dependent. As p increases the proportion 
of measured reflected signal increases and begins to dominate and the slopes 
of the plots change sign. As a result, the desirable type of scene from 
which to obtain data to process is one viewed through a thicker atmosphere 
with low ground reflectance. The L (0.) data for this type of scene has 

3 I 

greater curvature, increasing the amount of parameter independence, and reduc- 
ing the number of local minima in x space. The worst case is that for p 'v 0.5 
where SKYRAD could attempt to fix a complicated function to a set of data 
which is nearly a straight line. 
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Figures 20 and 21 are plots of the path radiance at the sensor, Lp, vs p 
for the models considered. 

An important characteristic to notice in Figures 18-21 is that for a 
given t and p value, the corresponding values of Ep and L p are uniquely de- 
termined. This fact will be utilized later in the procedure for selecting 
input choices for the a^ parameters and evaluating the results predicted by 
SKYRAD. 

Finally, a Isast squares straight line fit was made to l p (e.j) vs 
exp(-x sees..) determining a^ and a^ for all the wavelength optical depth 
cases considered and these were tabulated along with the computed Ep values. 
This permits estimation of Ep and L p for a given choice of X, r, and p. 

Fitting Procedure for Numerical Technique 

The model atmosphere radiative transfer calculations were used to provide 
unknowns for analysis by the non-linear regression method. An unknown was a 
data set of L £ ( 9^ ) values that was obtained from the model calculation output 
for an azimuth angle of 90°. The wavelength of the data was provided, but not 
the optical depth or reflectance value. This was done in an attempt to simu- 
late the realistic case of actual measured data. The goal was to utilize the 
numerical analysis method described to reproduce the initial t and p values 
of the model case and predict Ep and L . 

The experience gained analyzing the data of SYNRAD demonstrated that due 

2 

to the bumpy nature of the x hypersurface an arbitrary choice of input param- 
eters cannot be expected to converge to the correct result nor can a series of 
inputs be expected to converge to the same result. Therefore, a procedure was 
devised to judiciously determine a range of likely input parameters, to evalu- 
ate the subsequent outputs, refine the inputs and select a final result. 
Procedure: 

1) Select e 1 points for L g ( 0^ ) data and plot vs secG^ 

2) Compare plot of L s (e-) with curves from model calculations 

(a) From overall magnitude of L s (e^ ) estimate likely range of t, (a^ ) 

(b) From overall curvature of L g ( ) estimate likely range of p, (a 2 ) 

3) For range of t, p, from model calculation tables 
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Fig. 20. Plots of path radiance, L , versus reflectance, p, for 

the thin optical depth cases at the specified wavelengths. 
. Note the systematic increase in slope and magnitude. 


Plots of path radiance, Lp, versus reflectance, p, for 
the thick optical depth cases at the specified wavelengths 
Note the systematic increase in slope and magnitude. 





(a) Estimate range of Eq, (a j 

3 

(b) Estimate range of l p , (a^, a g ) 

4) Run SKYRAD for the selected range of input parameters 

5) Evaluate results and refine input parameters 

6) Run SKYRAD for the selected range of input parameters 

7) Select result 

Figure 22 illustrates the input and computed output L (©, ) values of this 

w I 

procedure for one of the unknowns considered. The selection of data points 
from the model calculations was not uniform, but concentrated points by in- 
terpolation in the ranged 15° < e.. < 45°, in order to emphasize the 
curvature in the input data. 

The initial range of parameter estimates wer§; 
x : 0.15 <a^ ■< 0.25 

p : 0.05 <a^ < 0*15 

E^: 0.085 <a g < 0.095 

L p : .070 <a^ < 0.075 

; .065 <a g < 0.070 

The criterion for evaluation of the SKYRAD parameter predictions were: 

2 

' 1. Is the reduced x value reasonable? 

2. Is the predicted value for x reasonable? 

> T RAYLEIGH ^ 

3. Is the value for E D reasonable? 

(SKYRAD tends to predict E g with a large uncertainty, but pathological 
cases E q > 0.35E Q or E D < 0.01E ) are still rejectable) 

4. Are the computed values of E g and L p consistent with the computed values 
for x, p? 

5. Langley plot (In (L (6.) - L (0. ) vs sec0.) the input data with the pre- 

^ » U 1 I 

dieted path radiance subtracted. Does the slope of this plot equal the 
predicted x value? • 

Approximately five sets of input parameters selected from the estimated 
ranges were run. The criteria above, were used to evaluate the resulting 
SKYRAD predictions and define three refined sets of parameter inputs. The 
outputs were again evaluated and the input parameters narrowed to those of 
the final result tabulated here and plotted in Figure 22. 
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Fig* 22. Plot of SKYRAD predicted fit to input model calculation 
unknown. Solid line represents radiative transfer 
calculation model data versus nadir angle. Dashed line 
is curve predicted by SKYRAD. 


FINAL RESULTANT PARAMETERS FOR MODEL DATA FIT 



T 

P 

e d 


' S 

Input 

.2 

.1 

.1 

.06 

.06 

Output 

.158 
+ .024 

.104 

+.023 

.124 
+ .153 

.061 
+ .004 

.059 
+ .005 

Model 

0.185 

0.1 

0.095 

0.062 

0.060 


Conclusion 

The non-linear regression numerical technique described is a viable 
method for determining the path radiance and ground reflectance in remote 
sensed scenes. The technique is not without inherent difficulty, most 
notably the problem of local minima. However, judicious application of 
the outlined analysis procedure and the careful evaluation of results 
yields parameters that are close to the correct values. 
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6.0 CONCLUSIONS & RECOMMENDATIONS 


The following conclusions are drawn from this study: 

1) The magnitude of the pointing angle from nadir should be restricted 
for atmospheric studies to no more than 50°. The reason is that no 
gain in accuracy is achieved in determining atmospheric parameters 
using large pointing angles, while in practice, the larger ground 
projected IFOV and the greater path radiance crosstalk from surround- 
ing ground areas may adversely affect the accuracy. In addition, 
engineering and multi -experiment requirements militate against the 
use of very large pointing angles. 

2) An increase in the atmospheric extinction coefficient, t, of 0.1 at 
0.65 ym caused an average decrease in classification accuracy of 3 %, 
Such a. change in x is large across a Landsat scene in a purely rural 
area. Thus atmospheric correction may not be an important task for a 
pointable imager for a rural scene. However, we should emphasize that 
we did not explore the case of low reflectances in bands 4 and 5 asso- 
ciated with very dense green vegetation. Also, this conclusion doss 
not apply to the important case of classification in the transitional 
zone between urban and rural areas. 

3) There is a significant asymmetry of classification accuracy between 
training and test data which indicates that, in general, it is prefer- 
able to train under hazy conditions and test under clear conditions 
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than vice versa. In addition, when training under clear and testing 
under hazy conditions, low reflectance features. Reversing the pro- 
cedure caused high reflectance features to show a greater error than 
low reflectance features. 

4) The two-wavelength approach to determining the extinction coefficient 
does not work under the conditions explored here. 

5) The non-linear least squares regression method proved to be success- 
ful in determining the extinction coefficient, the ground reflectance 
and the path radiance. However, care had to be exercised in prudently 
choosing the starting input parameters and in analyzing the output 
data. 

6) The analysis performed thus far does not indicate the need to extend 
the wavelength range beyond that required by other applications. 
Probably four bands in the range 0.45 - 0.75 pm should be sufficient to 
determine the spectral characteristics of the measured atmospheric 
parameters and to provide a check on the consistency of the results. 

7) Although the relative radiometric accuracy question was not addressed 
in this study, there are indications that 8 bit quantization may be 
marginal for techniques involving subtraction of radiance values, ' 
e.g., the two wavelength_ method. If the option of 9 or 10 bit quanti- 
zation exists* this question should be studied in detail. In any case 
future study should address the effect of quantization noise on class- 
ification accuracy. 

8) Although the correction for atmospheric effects may not often be 
beneficial for improving classification accuracy across purely rural 
scenes the opportunity to obtain scene reflectance values should be of 
great importance both for feature identification purposes, and BRDF and 
polarization studies . 

9) The results of this preliminary study are promising and indicate the 
feasibility of using a pointable imaging system (a) for the determina- 
tion of the atmospheric parameters required to improve classification 
accuracies in urban--rural transitional zones and to apply in studies 
of BRDF and polarization effects and (b) for the determination of the 
spectral reflectances of ground features. 
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The following recommendations are made based on the results of this 

preliminary 6-month study. (Of the six months, approximately the first 

three were devoted to a literature study published in another report). 

Future work should consider: 

1) Different combinations of solar zenith and azimuth angles. 

2) The separation of very similar scene classes (e.g. vegetation types) 
and classes having low reflectance values, 

3) The effect of detector noise on classification accuracy and reflectance 
estimation. 

4) The effect of the number of quantization steps on classification ac- 
curacy and reflectance estimation. 

5) The effect of an inhomogeneous atmosphere and a non-uniform ground 
scene, 

6) The effect of changing probability thresholding on classification 
accuracy. 

7) Typical natural variability differences between training and test areas. 

8) The very important question of whether very rapid, atmospherically in- 
duced irradiance fluctuations render greater than 6 bit quantization 
meaningless. 
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APPENDIX A 


A-l 


1 PROGRAM SKYRAD (INPUT, OUTPUT, TaPE 5,TAPE10, TAPEl 1) 

THIS ROUTINE UTILIZES THE NON-LINEAR REGRESSION SUBROUTINE CURFIT 
TO FIT A SET OF DATA TO A SPECIFIED FUNCTION WHICH IS NON-LINEAR 
IN ITS PARAMETERS, 

THE FUNCTION ANALYZED IS DEPENDENT ON FIVE PARAMETERS AND IS AN 
EXPRESSION FDR THE TOTAL RADIANCE RECIEVED BY A POINTABLE REMOTE 
SENSING SATELLITE-BASED PAOIOMET£R, (SEE SUBROUTINE FUNCT ) 

THE SOURCE OF SUBROUTINES CURFIT, FOERIV, AND SCHISQ IS 
BEV INGT0N,P,R , , DATA REDUCTION AND ERROR ANALYSIS FOR THE 
PHYSICAL SCIENCES, MCGRAW-HILL, 1969, 

description or parameters 
x j array nr independent variable csecthnj 

V » ARRAY OF DEPENDENT VARIABLE (RECEIVED TOTAL RADIANCE) 

SI KM AY ! ARRAY OF STANDARD DEVIATIONS FOR Y 
A t ARRAY OF fit PARAMETERS for FUNCTION 

OF.I.TAA l ARRAY OF INCREMENT VALUES FOR A 
SIGMAA I ARRAY OF STANDARD DEVIATIONS FOR A 
YFTT I ARRAY OF CALCULATED VALUES Of V 
NPTS t NUMBER OF PAIRS OF DATA POINTS 
NTEPHS t NUMBER DF FIT PARAMETERS 

MODE I DETERMINES METHOD OF STATISTICAL FIT (SEE CURFIT) 

CHI SOR » REDUCED CHI SQUARE VALUE FOR FIT 

Eo : EXTRATERRESTIAL solar IRRADIANCE VALUE FOR SPECTRAL band 
INPUT (UNFORMATTED READS) 

TAPES l CONTAINS X,Y,RIGMAY OUPUT OF SYNRAO OR REAL DATA 
T APE1P J NPT5,NTERMS,HODE,FLAMOA,A,DELTAA,EB 

CARDS t NFLG, IF NFLG * U PROGRAM READS T APE ID INPUT FROM CARD 

OUTPUT (UnFORMaTTFD WRITES! 

TApEH 1 NPTS,NTERMS,MOOE,FLAHDA, A,DELTAA 
PRINT l RESULTS OF EACH ITERATION THROUGH CURFIT 

COMMENTS 

1.) X,Y,SIGMAY, NPTS, MTERMS, MODE, FLAMDA, A, DELTA*, E0 MUST BE INPUT 
?.,) A,SIGHAA, YFIT,CH!SOR is OUTPUT 
3,) FLAMDA IS SET TO 0,001 INITIALLY 

«.) PROGRAM ITERATES THROUGH CURFIT UNTIL TIME IS CONSUMED 

OR DIFFERENCE OF CHISQR VALUES BETWEEN ITERATIONS IS < 0,001, 

IF CHISQR CONDITION IS NOT MET FINAL FIT RESULTS ARE STORED IN 
T APE l 1 TO BE USED FOR NEXT RUN OF SKYRAD, 

S3 PROPER COMMON BLOCK SPECIFICATION IS NECESSARY 
C h,) BEVIWGTON SUBROUTINES CURFIT, FOERIV, SCHISO ARE NEEDED 
C 

C THF. ROUTINE Was CODED IN FORTRAN for A CDC CYBEP 175, 

V* r BY S.J. MARTInEK, OPTICAL SCIENCES CENTER, OCTOBER 1979 

c 

•COMMON /ARRAY \ / X(l6), A(5), DELTA A (5) , 5IGHAA(5) 

COMMON /ARRAY?/ Y ( 1 6 ) , YFTT (16) , SIGHAY (16) 

CnMMON t / CHISQR, FLAHOA, E0 /NUM/ NPTS, NTERMS, MOOE 
55 DATA SIGMAA /5*0,0/ 

C 

C INPUT DATA 
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PROGRAM SKYPAO 70/7A 0PT«2 


60 


6 * 


7 t' 


75 


hr 


e r > 


90 


9 b 


100 


JPf 


nr 


READC5) X,Y,SXGMAY 
READ* , NFUG 
IF(NFLG) 2,2,3 

2 R£AD(10) NPYS,NTERM5,HOOE,FUAHOA,A,OELTAA,EB 
GO TO 0 

3 READ*, NPTS,NTERMS,MODE,FLAHDA,A,D£UAA,E0 
fl CONTINUE 

C 

C CALL CURFJT, PRINT RESULTS, AND ITERATE 
K • C 3 CHXSQR • e, 

PRINT SB 1 PRINT 55, K ,CHIS0R ,FLAMDA 
DO 5 X»1 , NTERMS 

5 PRINT 60, I.A(I),SIGHAa(I),DELTAA(I) 

CHITE5T * 1,0 
DO IS K* 1 , 35 

IF CMOO tK , fl) ,EO,0) PRINT «9 
CALL CURFIT 

PRINT 56, K,CNI3QR,FLAMDA 
DO IP 1*1, NTERMS 

10 PRINT 61, I , A (T ) , SIGMA A (I) 

IP CCHISQR , EQ , 0 , 0 5 GO TO 20 
C * ABSCCMISQR-CHITESTJ/CHITEST 
IF CS » 1 7 ,0,0101 SO TO 20 
CHITEST » CHISOR 
15 CONTINUE 
20 CONTINUE 

IF(K.E0,36) K ■ K-l 

C 

C OUTPUT RESULTS 

WRITE Cl 1 3 NPT3, NTERMS, MODE, FLAMDA, A, OELTAA^EB 
PRINT 65, K 
00 25 I*i,NPTS 

AE * Y(l)-YFITCI) 

PE * AE / Y ( I ) * 1 00 , 

Print 70, I,X£I),Y(!),YFXT(X),AE,P£ 

25 CONTINUE 
C 

C COMPUTE S PRINT PATH RaOIANCE FOR SECCTHN) ■ 1, 

PR AO * A ( a) « A ( 5 ) *EXP C»a Cl)) 

PRINT 75, PRAD 

c 

C FORMAT STATEMENTS 
fl9 FORMAT (IH1///) 

5e FORMAT (1H1//30X, "CURFIT ANALYSIS OF RADIANCE EXPRESSION V/13X 

2 ,"AJ ■ TaU > A? » RHO t A3 * EDIFF f A A , A5 * EXPANSION CONSTANTS 

3 " FOR PATH RADIANCE"/) 

55 FORMAT C1HB/2BX, "FIT ITER AT ION" , 13 , 5X , "CHISQR *",F6,3,5X, 

2 "FLAMOA 1PEA,1//56X, "A", 1 3 X , "DELA"/) 

56 FORMAT Cl H0/2EX, "FIT ITER ATION" , 13 , 5X , "CHI 50R *",F6,3,5X, 

? "FLAHOA »", 1PES. 1//51X, "A") 

60 FORMAT ( 1 H ,3BX,I5,3X,1PE11|4," E 1 1 AX ,»(", Ell ,fl,")" ) 

61 F0RMATC1H , 3BX, I5,3X , 1PE1 1 ,4," */.",Ell,A) 

65 FORMATC1HI///20X, "ERRORS IN RESULTS FOR ITERATION", I3//32X, 

? "X", 11X,"Y",SX,"YFTT",5X,"ABS ERR" , 6X , "PERCNT" /) 

70 FORMAT f 1 H ,2PX,I3,F9,?,AC1PE12.35) 

75 FOPMATUH0/15X,*»THE UPWELLING PATH RADIANCE IS",1PE12,3, 

2 " FOP A NAOJR ancle of ZERO") 
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subroutine 

3 

5 

10 

15 

20 

25 

3ft 

35 

5ft 
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CUHFIT 74/70 0PT«2 


SUBROUTINE CURFlT 

PURPOSE 1 MAKE A LEA3T»S0UaRES FIT TO A NON*LINEAR FUNCTION 
VIA THE MAROUUaROT APPROACH (SEE BEVlNGTON P,232) 

parameter description 

t ARRAY OF DATA POINTS FOR THE INDEPENDENT VARIABLE 

i array of data points for the dependent variable 
I ARRAY OF STANDARD DEVIATIONS for DEPENDENT varaible 
I NUMBER OF pairs of DATA P0INT5 
I NUMBER of FIT PARAMETERS 

t determines method of weighting least*squares fit 

♦ I (INSTRUMENTAL) HEIGHT (II ■ l , /SIGMAY (I ) **2 
0 (NO WEIGHTING) WF.IGHTCJ) * 1, 

*1 (STATISTICAL) HEIGHT (I) » 1,/Y(I) 
t ARRAY OF FIT PARAMETERS 

t ARRAY of increments for parameters a 
; array -of standard DEVIATIONS for PARAMETERS a 
1 PROPORTION OF GRADIENT SEARCH INCLUDED 
S ARRAY OF CALCULATED VALUES OF Y 
I REDUCED CHI SQUARE FOR FIT 
C 

C SUBROUTINES AND FUNCTION SUBROUTINES REQUIRED 
C FUNCT I EVALUATES THE FITTING FUNCTION 
C 5CHIS0 t EVALUATES REDUCED CHISOR FOR FIT TO DATA 

C FOERIV J EVALUATES THE DERIVATIVE OF THE FITTING FUNCTION 

C FOR ALL X VALUES FOR EACH PARAMETER 

C MXINV t INVERTS THE MATRIX 
C 

C COMMENTS 

c common block specifications APE CRUCIAL FOR parameters 
C DIMENSION STATEMNT VALID FOR NTERMS UPTO 10, NPT5 • U 
C SET FLAMDA ■ pi.eci AT BEGINNING OF SEARCH 

C NOTE 1 BEVJNGTON ROUTINE HAS BEEN ALTERED TO ONLY LOOK FOR POSITIVE P 
C 

COMMON /ARRAY1/ X(lfe), A(5), OELT A A (5) , SIGMAA(5) 

COMMON /ARRAY2/ V(i6), Y F I T (16) , SIGMAY (Ifc) 

COMMON / / CHISOR, FLAHOA, E0 /NUM/ NPTS, NTERMS, MODE 
COMMON /ARRAYS/ DERIV (10, lfe) , ARRAY (10, 10) /!/ NORDER 
DIMENSION HEIGHT(100), AUPMA (10, 10) , BETA (10), B(10), ADUM (18) 

C 

NORDER ■ NTFRMS 
HFPEE • NPTS^NTERMS 
IFCNFREE) 13,13,20 
13 CHISQR ■ 0, 

GO TO 110 
C 

C EVALUATE WEIGHTS 
2P DO 30 I ■ 1 , NPTS 



IF (MODE) 

IZi 

'27,29 

22 

IF(YCI)) 

H5, 

,27,23 

23 

WEIGHT (I) 

a 

l./YCI) 


GO TO 30 



25 

WEIGHT (I) 

« 

i./(*Y(in 


GO TO 30 



27 

WEIGHT ( I ) 

a 

i. 


X 
Y 

SIGMAY 
NPTS 
L NTERMS 
C MODE 
C 
C 

c 

C A 

c deltaa 

C SIGKAA 
C FLAMDA 
C YFIT 
C CHI5DR 
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SUBROUTINE CURFXT 70/70 OPt«g 


61 ' 


66 


7P 


76 


6 0 


85 


9r 


9** 


10P 


1 PIN 


UP 


CO TO 30 

29 WeiGHY(I) ■ i,/8XCH*y(U«*2 

30 CONTINUE 
C 

C EVALUATE alpha and beta MATRICES 
00 34 J® 1 » NTERMS 
BETACJJ • 0, 

00 30 K«i,J 

3« ALPHA (J , K) » 0, 

CALL FDERIV 
CALL FUNCT 
DO 50 It! f NPT3 
DO 46 J»1 , NTERM3 

BETA (.T) ■ BETA (J) ♦ HEIGHT (I) * CY (!) -YFIT (I )) *DERIV (J , I) 

00 06 KiUJ 

06 ALPHA(J,K) • ALPHA ( J* K) ♦ WEIGHT Cl) *OERIV (J ,!) *OERIV (K , I) 

50 CONTINUE 

DO 53 J® 1 , NTERM5 
DO 53 K» U J 

53 ALPHA f K , J) • ALPHACJiK) 

C 

C EVALUATE chi SQUARE AT STARTING POINT 
CALL 3CHIS0 
CHISDi ■ CHISOR 
C 

C INVERT MOnTFlEO CURVATURE MATRIX TO FINO NEW PARAMETERS 
KK ■ B S NN « 0 
71 M I KK4-1 

DO 7 0 jm,NTERMS 
PO 73 K»1,NTERMS 

73 ARRAY ( J , K 3 ■ ALPHA(J,K)/S0RT(ALPHA(J,J)#ALPHA(K>K)7 
7 0 ARRAY (J » J3 ■ 1,*FLAM0A 
CALL M X T N V 
DO M J ■ 1 , NTERHf 
B(J3 * ADUM(J) a A ( J 3 
DO 80 K«irNTERHS 

60 6 (JJ ■ B(J)+BETA (K)*ARRAY ( J , K) /SORT (ALPHA (J , J )* ALPHA (K , K ) ) 

DO 85 J* U NTERH5 
05 A ( J ) • 6 C JD 
C 

C IF CHT SQUARE INCREASED t INCREASE FLAMOA AND TRY AGAIN 
CALL FUNCT 
C ALL 5CH I SO 

'PRINT 150, KK,FLAMDA,CHISOi,CHISQR 

150 FP>RHATUHB,UX,»LAMDA ITER A T I ON* , 1 3 , 3X , "FL AMD A ■ " , 1PE6 , 1, 3X , 

2 "CHI SO l »",0PFB, 3, 3X , "CHISOR ■",F6,33 

C 

C CHECK IF PARANETER5 HAVE GONE NEGATIVE AND TRY TO FORCE POSITIVE 
NN • NW+1 

XX » AMlNUA(l),AC2),AC3),AC<n,AC5)) 

IF (XX ) 66,86,90 
66 PRINT 151, NN 

15 1 FORMAT (1H0,23X, "NEGATIVE PARAMETER ITERATION" , 133 
IFCNN.LT. 73 GO To 95 

PRINT 152 

152 FORHAT(1H0,7X, "FIT TERMINATION t CANNOT AVOID OPTIMIZING A" 

2 " NEGATIVE PARAMETER VALUE") 
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A-S 


SUBROUTINE 

US 

12 " 

12 '< 


0 0 «L PflQE |s 
°F POOR QUALITY 


CURFJIT 7a/?fl 0FT«8 


CHISOR > 0,0 
GO TO 101 

90 IP (CHISQUCHISOR) 95.101*101 

95 FL*MOA • i&,*Fl*HOA 
00 96 Jul.NTERHS 

96 A ( J) • ADUMtJ) 

GO TO 71 

EVALUATE PARAMETERS and UNCERTAINTIES 
101 on 103 J*l . NTERMS 
i 03 SIGMAA(J) » .SGRT(ARRAT(J, J)/ALPHA(J, ju 
FLAMDA » FLAMDA/10, 

110 RETURN 
£NQ 
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Si*y»ntIT JNE FUNCT 70/70 0PT»2 


1 




\5 


2 '* 


25 


3^ 


SUBROUTINE FUNCT 
C 

C THIS ROUTINE COMPUTES YFIT , FOR ALL X, FOR A GIVEN CALL 
C 

C YFIT (I) • A* ♦ C(E0*EXP(«At*SEC7HZ) ♦ A3)*A2/P3 * A 55 »EXP C-A 1 #X (I ) ) 

C 

r PflWQT&ATLlTA t 

c kVa solar zenith angle of osdeg is assumed 

C 2,1 E55 ■ EP*CDS(THZ) 

C 

COMMON /ARRAY 1 / X ( 16) « AC5)» DELTAA (5) » SIGKAAC5) 

COMMON /ARRAy2/ YCt<0, YFTTCI6), SIGMA Y C 1 6 J 
COMMON / / CHISOR, FLAMOA, E0 /NUM/ NPTS? NTERMS, MODE 
C 

PI * ACOSC-1.) S RPD • PI/180, » TH2 » 05, 

A l ■ A ( ] ) S A2 ■ A f,2) S A3 * ACS) S A4 » A(0) S A5 ■ *C5) 
C 

n compute terms of expression 

ET « EXP(-Al) 

RMZ • l f /caS(TMZ*RPO) 

PI s EBs£T**RMZ/RMZ 
p? * CP1*A3)*A2/PI - *5 
C 

C FILL YPJT ARRAY 

DO § !■ I # NPTS 

YFIT Cl) ■ AO t P2*ETt*XCI) 

5 CONTINUE 

return 

END 
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A-7 


1 


*5 


1# 


1? 


pi* 


?s 


3 .A 


Vi 


subroutine fderiv 

C PURPOSE l EVALUATE DERIVATIVES OP FUNCTION FOR LEAST-80R AURES SEARCH 
C AT EACH POINT X WRT EACH FIT PARAMETER 
C 

C FOR ARBITRARY FUNCTION GIVEN BY F'JNCT 

C 

C SEE BEVINGTON, Pgfl2 
C 

c description of parameters 

C X t ARRAY OF OATa POINTS FOR INDEPENDENT VARIABLE 

C I t INDEX OF DATA POINTS 

C A : apray OF parameters 

C DELTAA r ARRAY OF PARAMETER INCREMENTS 

C NTERMS j NUMBER OF PARAMETERS 
C OERIV | DERIVATIVES OF FUNCTION 
C 

COMMON /ARRAY 1 / X C 1 6 3 > A C53 , DELTAA (S3 » SICMaA C 5 > 

COMMON /ARRAYS/ Y(163» YFITClbJ « SIGMAY (lfc) 

COMMON / / CHJSOR, FLAMOA, EO /NUM/ NPTS, NTERMS, MODE 
COMMON /ARRAYS/ DERIVdfi, i<>3 , ARR A Y (J Pi , 1 03 / 1 / NORDER 
DIMENSION YH323 
C 

DO 15 Ja l , NTERMS 

A J a A C J3 S DEL ■ CELT AA ( J) 

A C J3 a AJ * DEL 
call PUNCT 
OO 5 I ■ 1 ,NPTS 
5 Y1CI3 ’» YFIT(I) 

AtJ3 « AJ-OEL 
CALL FUNCT 
DO ID I«1,NPT5 

IP OERIV CJ » j3 « t Y I CI3-YFITCI))/f2,»OEL) 

A C J3 ■ AJ 
15 CONTINUE 
RETURN 
END 
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SUBROUTINE 
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1* 

1 S 
2"' 

7C, 

3S 
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SCHISa 7U/7U OPT*g 


SUBROUTINE SCHISO 

EVALUATES The REDUCED CHI SQUARE FOR FIT TO DATA 
CHISOR ■ SUM ( (Y-YFIT) **2/SIGMA**2) /NFREE 
SEE REVINGTON, PJ93 

USAGE t PROPER COMMON BLOCK SPECIFICATION 

description OF parameters I 

Y - ARRAY OF DATA POINTS 

5IGMAY • ARRAY OF STANDARD deviations FOR DATA POINTS 

NPTS - NUMBER OF DATA POINTS 

NFREE • NUMBER OF DEGREES OF FREEDOM 

MODE • DETERMINES METHOD OF WEIGHTING LEAST^SOliARES FIT 

ft (INSTRUMENTAL! WEIGHT(I) ■ l , /SIGMA Y ( ") **2 
0 (NO WEIGHTING) wEIGHT(I) ■ 1, 

-1 (STATISTICAL) WEIGHT (I ) * t./Y(T) 

YF1T - ARRAY OF CALCULATED VALUES OF Y 
CHISOR - REDUCED CHI SQUARE FOR FIT 


COMMON /(ARRAY l / X(16), A(5), DELTAA(5 )i SIDMAA(5) 
COMMON /ARRAY2/ Y (16) • YFIT(lb), SIGMAY ( 1 b) 

COMMON / / CHISOR. FLAMDA, ED /NUM/ NPTS, NTERMS , MODE 


CHISQ • 0, 

NFREE ■ NPTS-NTERMS 
IF(nFUEE) 13,13,20 
13 CHISOR » 0. 

GO TO «H 

ACCUMULATE CHI SQUARE 


20 DO 30 I « 1 , NPTS 
IF (MODE) 22,27 ,29 

22 I F ( Y ( I ) ) 25,27,23' 

23 WEIGHT ■ 1 , / V ( I) 

GO TO SO 

25 WEIGHT ■ 1 / ("Y (I ) ) 

GO TO 3D 
27 WEIGHT * 1, 

GD Tfj 30 

29 HEIGHT ■ 1,/SIGMAV(IJ**2 

30 CHISQ • CHISQ ♦ WEIGHT* (Y CI)*YFIT (I) ) *»2 

C DIVIDE BY NUMBER OF DEGREES OF FREEDOM 
C 

FREE * NFREE 
CHISQR ■ CHISO/FREE 
A0 RETURN 
END 
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SUBROUTINE MXJ MV 7il/7 A 0PT»2 


5 


5 


if* 


IS 


if 


25 


3P 


35 


ap 


uf 


5f. 


SUBROUTINE mxinv 

c 

C THE OBJECTIVE OF THIS ROUTINE IS TD COMPUTE THE INVERSE OF HARR I X 
C ARRAY, THE METHOD USED IS GAUSS. JORDAN , NO ROW SUBSTITUTIONS, 

C THE MATRIX ARRAV IS ASSUMED TO BE NON-SINGULAR AND THE DIAGONAL 
C ELEMENTS ARE USED AS PIVOTS, 

c note hi the original matrix array is destroyed, the inverse is 

c returned in MATRIX ARRAY 

c 

COMMON /ARRAYS/ DERI v (l<M&n ARRAY C 1 (S « 105 / 1 / NORDER 

real invcu, h») 

c 

C PUT IDENTITY matrix IN INVERSE MATRIX ARRAY 
N • NORDER 
DO ?p I » 1 * N 
DO 25 J ■ 1 1 N 
INVCI.J) • 0, 

25 CONTINUE 

iNV(i,n ■ i, 

20 CONTINUE 
C 

e find inverse using gauss. Jordan method 

DO 50 I»1|N 
C 

c select pivot element 

ALPHA ■ ARRAY (IiH 
C 

C DIVIDE PIVOY ROW BY PIVOT ELEMENT 
DO 30 J ■ 1 , N 

APR AY ( J , J ) s ARRAY (I, J) /ALPHA 
INVCHJJ * INV (I # J5 /ALPHA 
30 CONTINUE 
C 

C FOR DTHFR rows, SUBTRACT MULTIPLIER TIMES PIVOT ROW FROM ROW 
DO U0 K * 1 , N 

IF CK.EO.I5 GO TO A0 
BETA * ARRAYOHI) 

DO 35 L*WN 

APRAYCK,L3 * ARRAY CK ,L)-BETA«APRAT (IrL5 
I N V (K , L I » INVcK,L)-BETA»INVCI,L) 

35 CONTINUE 
A0 CONTINUE 
SB. CONTINUE 
C 

C FILL ARRAY WITH INVERSE VALUES 
00 *0 J » 1 , N 
DO 60 1*1 ,N 

ARRAY ( I , J) « INV(I,J3 
60 CONTINUE 
RETURN 
END 
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APPENDIX B 

SAMPLE SKY RAD OUTPUT 

CUSFJT ANALYSIS OF RADIANCE EXPRESSION 

A 1 8 TAu > A? ■ RHP f A3 » ED IFF » AA,A5 ■ EXPANSION CONSTANTS FOR PATH RADIANCE 


FIT ITERATION 0 CHI8QR * 0,000 FLAMDA a 1.0E-03 


l, AN l,t A ITERATION 


7 f 5000E-02 +/• 0, 
8.S000E-01 + /- 0, 
5.0030E-02 */• 0, 
7.50P0E-02 */- 0, 
7.5000E-02 */* 0, 


( l,0000E-02) 
( i , 0000E*02) 
C i,0000E-02) 
( l,0B00E*02) 
C 1.0000E-02) 


FLAMDA * 1 , 0E-03 CHISD1 *2038,036 CHISQR 


1,261 


FIT ITERATION 


CHISQR • 1,261 


FLAMDA « 1.0E-0A 


j.AHOA iteration 


1.1223E-01 */- 1.C793E-02 
2.5705E-P1 ♦ A ,2369E»03 
6.9344E-02 */• 1.1625E-02 
7,b956E-02 +/• 8,0684E-0fi 
7,3fl61E-02 */• 9,c?5l0E-0« 

FLAMOA e 1.0E-04 CHISQ1 ■ 1,261 


CHI SDR 


FIT ITERATION 


CHISQR ■ 


FLAMOA ■ l.BE-05 


1, a099E.pl ♦ /- 2.0926F-02 
2.5813E-01 +/• 1.A2BPE-02 
7.2170E-02 */• 3.T37AE-02 
7 , 763SE-02 ♦ /« 2.3B4SE-03 
7.3259E-02 +/• 3.0S60E-03 


L AMP A ITERATION 1 FLAMOA ■ 1.0E-0S 


CHJSQ1 


,7b« CHISOR * 


FTT ITERATION 3 


CHISQR * 


FLAMDA « 1 , 0E-06 


1.3761E-01 ♦ /- 6,3818E-02 
2 ,56946*01 + /• 4.B223E-02 
6,917 3E-02 ♦ 1.217BE-01 
7.6078E-02 + /- 7.4089E-03 
7.3476E-02 ♦ /- 9,99986-03 
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B-2 


LAHDA ITFRATIDN 1 
L AMD A ITERATION 2 


FLAMDA • 1.0E-06 
FLAMDA » i,0E*05 


CHX3Q1 ■ ,444 

CHISOI » ,444 


CHI30R * 1,636 

CHI80R • ,444 


FIT ITERATION 


CHI3QR • ,444 


FLAMDA » 1.0E-06 


1 

2 

3 

4 

5 


i‘,3420E»0l 
2,5588E»ol 
6 .6468E-02 
7 , 85 1 7E-0? 
7 ,3697E®02 


A 

*/- 6.106BE-02 
♦/- 4,601 1£»02 

♦ 1 ,2166E»01 
*/• 7.4740E-03 

♦ /- 9.9523E-03 


ERRORS IN RESULTS FOR ITERATION 4 



X 

V 

YF IT 

ABS ERR 

PERCNT 

1 

1,00 

6.063E-02 

b,04bE»02 

1,61 6E»04 

Zf 669E«0 1 

2 

1,00 

6.056E-02 

6.047E-02 

9.066E-05 

1.501E-01 

3 

1.01 

6.B52E-02 

6.049E-02 

3.300E-0S 

5.453E*02 

4 

1.P2 

6.051E-02 

6.052E-02 

•1.225E-05 

«2 , 025E»02 

S 

1,04 

6.052E-02 

6.056E-02 

-4.05AE-05 

•6.705E-02 

6 

1 ,06 

6.0S6E-0? 

6.B62E-02 

-5.803E-05 

»9, 582E*02 

7 

1,09 

6.063E-02 

6.069E-02 

•6 , 6OflE*05 

-U0B9EW01 

B 

1.53 

6.072E-02 

6.076E-02 

-6.156E-05 

•t,014E»01 

9 

1.15 

6.07BF-02 

6.0B3E-02 

*5,4l(,E.05 

•6,91 0E*02 

10 

1,16 

6.0B3E-02 

6.0B9E-02 

• 6 , 2 15E»05 

• 1 .B22F-01 

11 

1,21 

6.091E-02 

b,096E-02 

-4 .599E-05 

•7.551E-02 

12 

1,24 

6, 100E-02 

6, 103E-0? 

-2.620E-05 

-4.295E-02 

13 

1.27 

6, H0E-O2 

6.110E-02 

•3.372E-06 

•5.516E-03 

1« 

1 ,31 

6. 120E-02 

6, 1 19E-02 

1 , 160E-05 

1.927E-02 

15 

1.35 

6.133E-02 

6,128E*02 

4,850E*O)5 

7 , 908Ei»02 

16 

1.39 

6, 149E-02 

6.136E-02 

1 ,056E*04 

1 ,721E»01 

UPWEIL1NG 

PATH 

RADIANCE IS 

1 .406E-02 

FOR A NADIR 

angle of zero 
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SAMPLE OUTPUT OF RADIATIVE TRANSFER COMPUTATION 
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